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Abstract. This paper introduces a new weak Galerkin (WG) finite element method for second 
order elliptic equations on polytopal meshes. This method, called WG-FEM, is designed by using 
a discrete weak gradient operator applied to discontinuous piecewise polynomials on finite element 
partitions of arbitrary polytopes with certain shape regularity. The paper explains how the numerical 
schemes are designed and why they provide reliable numerical approximations for the underlying 
partial differential equations. In particular, optimal order error estimates are established for the 
corresponding WG-FEM approximations in both a discrete H 1 norm and the standard L 2 norm. 
Numerical results are presented to demonstrate the robustness, reliability, and accuracy of the WG- 
FEM. All the results are derived for finite element partitions with polytopes. Allowing the use of 
discontinuous approximating functions on arbitrary polytopal elements is a highly demanded feature 
for numerical algorithms in scientific computing. 
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1. Introduction. In this paper, we are concerned with a further and new devel- 
opment of weak Galerkin (WG) finite element methods for partial differential equa- 
tions. Our model problem is a second-order elliptic equation which seeks an unknown 
function u = u(x) satisfying 

(1.1) - V • (a(x,u, Vu)Vw) = f(x), in ft, 

where f2 is a polytopal domain in R (polygonal or polyhedral domain for d = 2, 3), 
V« denotes the gradient of the function u = u{x), and a = a(x, u, Vit) is a symmetric 
d x d matrix-valued function in Q. We shall assume that the differential operator is 
strictly elliptic in f2; that is, there exists a positive number A > such that 

(1.2) ^a(x,r],p)^ > A£*£, V£ e R d , 

for all x € ft, r\ S R,p £ R d . Here £ is understood as a column vector and £* is the 
transpose of £. We also assume that the differential operator has bounded coefficients; 
that is for some constant A we have 

(1.3) |a(x,i7,p)| < A, 

for all x E ft, r) e R, and p e R d . 
Introduce the following form 

(1.4) a(4>;u,v) := / a{x,<t>,V4>)Vu -\7vdx. 

Jo, 
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For simplicity, let the function / in (jl.ll) be locally integrable in SI. We shall consider 
solutions of (jl.ip with a non-homogeneous Dirichlet boundary condition 

(1.5) u — g, on 90, 

where g € H^{d£i) is a function defined on the boundary of SI. Here iJ 1 (SI) is 
the Sobolev space consisting of functions which, together with their gradients, are 
square integrable over SI. H^(dVl) is the trace of _ff 1 (Sl) on the boundary of SI. The 
corresponding weak form seeks u £ H 1 ^) such that u = g on <9S1 and 

(1.6) a(u;u,v) = F(v), Vu e flo(fi), 

where F(v) = J n fvdx. 

Galerkin finite element methods for (jl.6l) refer to numerical techniques that seek 
approximate solutions from a finite dimensional space Vh consisting of piecewise poly- 
nomials on a prescribed finite element partition Th ■ The method is called conforming 
if Vh is a subspace of i/ 1 (SI) . Conforming finite element methods are then formulated 
by solving Uh € Vh such that Uh = Ihg on 5S1 and 

(1.7) a(u h ;u h ,v) = F(v), Vv eV h nH%(p,), 

where Ihg is a certain approximation of the Dirichlet boundary value. When Vh is not 
a subspace of ff 1 (S7), the form a((f>;u, v) is no longer meaningful since the gradient 
operator is not well-defined for non-£f 1 functions in the classical sense. Nonconforming 
finite clement methods arrive when the gradients in a(</>; u, v) are taken locally on each 
element where the finite clement functions are polynomials. More precisely, the form 
a{4>;u, v) in nonconforming finite element methods is given element-by-element as 
follows 

(1.8) ah(4>',u, v) := / a(x, <f>, V<?!>)Vu • Vvdx. 

T&T h Jt 

When Vh is close to be conforming, the form dh(4>;u,v) shall be an acceptable ap- 
proximation to the original form a(4>; u, v). The key in the nonconforming method is 
to explore the maximum non-conformity of Vh when the approximate form a^(0; u, v) 
is required to be sufficiently close to the original form. 

A natural generalization of the nonconforming finite element method would occur 
when the following extended form of (|1.8j) is employed 




where V tt , is an approximation of V locally on each element. By viewing V w as 
a weakly defined gradient operator, the form a w ((j);u,v) would give a new class of 
numerical methods called weak Galerkin (WG) finite element methods. 

In general, weak Galerkin refers to finite element techniques for partial differential 
equations in which differential operators (e.g., gradient, divergence, curl, Laplacian) 
are approximated by weak forms as distributions. In |17j . a WG method was in- 
troduced and analyzed for second order elliptic equations based on a discrete weak 
gradient arising from local RT [TB] or BDM [5] elements. Due to the use of the RT 
and BDM elements, the WG finite element formulation of [T7] was limited to classical 
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finite element partitions of triangles (d — 2) or tetrahedra (d — 3). In [18] . a weak 
Galerkin finite element method was developed for the second order elliptic equation 
in the mixed form. The use of a stabilization for the flux variable in the mixed formu- 
lation is the key to the WG mixed finite element method of [18]. The resulting WG 
mixed finite element schemes turned out to be applicable for general finite element 
partitions consisting of shape regular polytopes (e.g., polygons in 2D and polyhedra 
in 3D), and the stabilization idea opened a new door for weak Galerkin methods. 

The goal of this paper is to apply the stabilization idea to the form a w (4>; u, v), and 
thus to develop a new weak Galerkin method for (|l.l[) - f|1.5[) in the primary variable 
u that shall admit general finite element partitions consisting of arbitrary polytopal 
elements. The resulting WG method will no longer be limited to RT and EDM 
elements in the computation of the discrete weak gradient V m . In practice, allowing 
arbitrary shape in finite element partition provides a convenient flexibility in both 
numerical approximation and mesh generation, especially in regions where the domain 
geometry is complex. Such a flexibility is also very much appreciated in adaptive mesh 
refinement methods. 

The main contribution of this paper is three fold: (1) the WG finite element 
method to be described in section [4] allows finite element partitions of arbitrary 
polytopes which are shape regular in the sense as defined in [IS], (2) the finite el- 
ement spaces constitute regular polynomial spaces on each clement/face which are 
computation-friendly, and (3) the WG finite element scheme retains the mass conser- 
vation property of the original system locally on each element. 

One close relative of the WG finite element method of this paper is the hybridiz- 
able discontinuous Galerkin (HDG) method [12 . In fact, it can be proved that our 
weak Galerkin method is identical to HDG method for the Poisson equation. How- 
ever, the WG method differs from HDG for the model problem (|l.ll) with nonconstant 
coefficient matrix a and more sophisticated problems. These two methods are funda- 
mentally different in concept and formulation. The key element of HDG is the flux 
variable, while the key element for WG is the gradient operator through weak deriva- 
tives. For either nonlinear or degenerate coefficient matrix a = a(x,u, Vit), the WG 
finite clement method has obvious advantage over HDG since Vit is approximated 
by S7 w u and there is no need to invert the matrix a in WG formulations. More im- 
portantly, the concept of weak derivatives makes WG a widely applicable numerical 
technique for a large variety of partial differential equations which we shall report in 
forthcoming papers. 

The paper is organized as follows. In section [2] we introduce some standard 
notations in Sobolev spaces. In section |3l we review the definition and approximation 
of the weak gradient operator. In section 01 we provide a detailed description for the 
new WG finite element scheme, including a discussion on the element shape regularity 
assumption. In section O we define some local projection operators and then derive 
some approximation properties which are useful in error analysis. In section [51 we 
show that the WG finite element method retains the mass conservation property of the 
original system locally on each element. In section [3 we show that the weak Galerkin 
finite clement scheme for the nonlinear problem has at least one solution. The solution 
existence is based on the Leray-Schauder fixed point theorem. In section [8) we shall 
establish an optimal order error estimate for the WG finite element approximation 
in a H 1 - equivalent discrete norm for the linear case of (jl.ip . We shall also derive 
an optimal order error estimate in the L 2 norm by using a duality argument as was 
commonly employed in the standard Galerkin finite element methods [5] . Finally 
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in section |H1 we present some numerical results which confirm the theory developed 
in earlier sections. 

2. Preliminaries and Notations. Let D be any domain in M. d ,d — 2,3. We 
use the standard definition for the Sobolev space H 3 (D) and their associated inner 
products (•, -) s ,d, norms || • \\ s .d, and seminorms | • \ Si d for any s > 0. For example, 
for any integer s > 0, the seminorm | • \ St ]j is given by 




with the usual notation 

d 

a = (ai, . . . ,a d ), \a\ = ai + . . . + ad, 3 a = JJ<9";. 

i=i 

The Sobolev norm || • \\ m ,D is given by 

\Hm,D= (X>IL^ ■ 

The space H (D) coincides with L 2 (D), for which the norm and the inner product 
are denoted by and (-,-)d, respectively. WhenD = Q, we shall drop the subscript 
D in the norm and inner product notation. 

The space H(div; D) is defined as the set of vector-valued functions on D which, 
together with their divergence, are square integrable; i.e., 

£T(div; D) = {v : v G [L 2 (D)] d , V • v e L 2 (D)} . 

The norm in iJ(div; D) is defined by 




3. Weak Gradient. The key in weak Galerkin methods is the use of discrete 
weak derivatives in the place of strong derivatives in the variational form for the un- 
derlying partial differential equations. For the model problem (jl.6l) . the gradient V 
is the principle differential operator involved in the variational formulation. Thus, 
it is critical to define and understand discrete weak gradients for the corresponding 
numerical methods. Following the idea originated in |17j . the discrete weak gradi- 
ent is given by approximating the weak gradient operator with piecewise polynomial 
functions; details are presented in the rest of this section. 

Let K be any polytopal domain with boundary dK. A weak function on the 
region K refers to a function v = {vo,vt} such that vo G L 2 (K) and Vb £ H?{dK). 
The first component vq can be understood as the value of v in K, and the second 
component Vb represents v on the boundary of K. Note that Vb may not necessarily 
be related to the trace of vq on dK should a trace be well-defined. Denote by W(K) 
the space of weak functions on K] i.e., 



(3.1) W(K) := {v = {v ,v b } : V G L 2 (K), v b G H$(dK)}. 
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The weak gradient operator, as was introduced in |17j . is denned as follows. 

Definition 3.1. The dual of L 2 (K) can be identified with itself by using the 
standard L 2 inner product as the action of linear functionals. With a similar inter- 
pretation, for any v G W(K), the weak gradient of v is defined as a linear functional 
S7 w v in the dual space of H{div,K) whose action on each q G H{div,K) is given by 

(3.2) (V m j), q) K ■= -(vo, V • q) K + (vb, q ■ n)g K , 

where n is the outward normal direction to dK, (vq, V • q)x — Jr- u o(V ■ q)dK is the 
action of vq on V ■ q, and (vb,q ■ n)dK is the action of q ■ n on Vb G H^(dK). 

The Sobolev space H 1 ^) can be embedded into the space W{K) by an inclusion 
map iw ■ H X {K) — s> W(K) defined as follows 

iw(<t>) = {4>\ K ,<P\dKh 4> e H l (K). 

With the help of the inclusion map iw, the Sobolev space H 1 ^) can be viewed as a 
subspace of W(K) by identifying each (f> G if 1 (If) with iw {<!>)■ Analogously, a weak 
function v = {vo,vi,} G W(K) is said to be in iJ 1 (if) if it can be identified with 
a function <f> G if 1 (if) through the above inclusion map. It is not hard to see that 
the weak gradient is identical with the strong gradient (i.e., V w v — Vw) for smooth 
functions v G ^(K). 

Recall that the discrete weak gradient operator was defined by approximating 
Vu, in a polynomial subspace of the dual of H(div, K). More precisely, for any non- 
negative integer r > 0, denote by P r {K) the set of polynomials on K with degree no 
more than r. The discrete weak gradient operator, denoted by V w ,r,K, is defined as 
the unique polynomial {\7 w ,r,Kv) G [P r {K)] d satisfying the following equation 

(3.3) (V w ^ K v,q) K = -(«o,V -q) K + {v b ,q-n) 9K , Vq G [P r (K)] d . 

The discrete weak gradient operator, namely ^ w ,r,K as defined in (|3.3[) . was first 
introduced in [17j where two examples of the polynomial subspace [P r (K)] d were 
thoroughly discussed and employed for the second order elliptic problem (|l.lj) - (|1.5p . 
The two examples make use of the Raviat-Thomas [16] and Brezzi-Douglas-Marini 
[5] elements developed in the classical mixed finite element method. As a result, the 
corresponding WG finite element method of [17] is closely related to the mixed finite 
element method. In this paper, we shall allow a greater flexibility in the definition 
and computation of the discrete weak gradient operator V Wt r,K by using the usual 
polynomial space [P r {K)] d . This will result in a new class of WG finite element 
schemes with remarkable properties to be detailed in forth coming sections. 

4. Weak Galerkin Finite Element Schemes. In finite element methods, 
mesh generation is a crucial first step in the algorithm design. For the usual finite 
element methods [111 [7], the meshes are mostly required to be simplices: triangles or 
quadrilaterals in two dimensions and tetrahedra or hexahedra in three dimensions, or 
their variations known as isoparametric elements. Our new weak Galerkin finite ele- 
ment method is designed to be sufficiently flexible so that general meshes of polytopes 
(e.g., polygons in 2D and polyhedra in 3D) are allowed. For simplicity, we shall refer 
the elements as polygons or polyhedra in the rest of the paper. 
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4.1. Domain Partition. Let Th be a partition of the domain il consisting of 
polygons in two dimensions or polyhedra in three dimensions satisfying a set of con- 
ditions to be specified. Denote by Eh the set of all edges or flat faces in Th, and let 
£ ° = £h\dfl be the set of all interior edges or flat faces. For every element T G Th, we 
denote by \T\ the area or volume of T and by Kt its diameter. Similarly, we denote 
by |e| the length or area of e and by h e the diameter of edge or flat face e G Eh- We 
also set as usual the mesh size of Th by 

h = max Kt. 

Ten 

All the elements of Th are assumed to be closed and simply connected polygons or 
polyhedra. We need some shape regularity for the partition Th described as follows 
(see [T5] for more details). 

Al: Assume that there exist two positive constants g v and g e such that for every 
element T £ Th we have 

(4.1) g v h d T < \T\, q^- 1 < \e\ 

for all edges or flat faces e of T. 
A2: Assume that there exists a positive constant n such that for every element 
T G Th we have 

(4.2) nh T < h e 

for all edges or flat faces e of T. 

A3: Assume that the mesh edges or faces are flat. We further assume that for every 
T G Th, and for every edge/face e G <9T, there exists a pyramid P(e,T, A e ) 
contained in T such that its base is identical with e, its apex is A e G T, and 
its height is proportional to hr with a proportionality constant a e bounded 
from below by a fixed positive number a*. In other words, the height of 
the pyramid is given by a e hx such that a e > a* > 0. The pyramid is also 
assumed to stand up above the base e in the sense that the angle between 
the vector x e — A e , for any x e G e, and the outward normal direction of e is 
strictly acute by falling into an interval [0,0q] with Oq < J . 

A4: Assume that each T G Th has a circumscribed simplex S(T) that is shape 
regular and has a diameter hs(T) proportional to the diameter of T; i.e., 
h>S{T) < 7*/it with a constant 7* independent of T. Furthermore, assume 
that each circumscribed simplex S{T) intersects with only a fixed and small 
number of such simplices for all other elements T Glu- 

Figure 14.11 is a depiction of a shape-regular polygonal element in 2D. As to the 
property A3, for edge e = AF, the corresponding pyramid is given by the triangle 
AA e F which is of a similar size as the polygonal element. n e is the outward normal 
direction to the edge e. The angle between the two vectors n e and A e x e is strictly 
acute for any x e G e. 

4.2. WG Finite Element Algorithms. Let Th be a finite element partition 
that is shape regular; namely, satisfying the properties A1-A4. On each element 
T G Th, we have a space of weak functions W(T) defined as in Section [3J Denote by 
V the weak function space on Th given by 



(4.3) 



V := {v = {v , v b } : {v , v b }\ T G W(T),T G %}, 
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Fig. 4.1. Depiction of a shape-regular polygonal element ABCDEFA. 

where {vo,vi,}\t '■— {(vo)|t, (vb)\dT} is the restriction of v on the element T. 

For any given integer k > 1, let Wk(T) be the discrete weak function space 
consisting of polynomials of degree k in T and piecewise polynomials of degree k on 
<9T; i.e., 

(4.4) W k {T) {v = {v a ,v b } : v \ T G P fc (T), v 6 | e G P fe (e), e £ 5T}. 
Furthermore, let Vh be the weak Galerkin finite element space defined as follows 

(4.5) V h := = {v ,v b } : {« ,« 6 }|t G W k {T), T G 77J 
and 

(4.6) := {v : v b = on <9ft}. 

Denote by V^fc-i the discrete weak gradient operator on the finite element space 
Vh computed by using (13.31) on each element T; i.e., 

(V w ,k-iv)\T = V^^tHt), Vi> G V h - 

For simplicity of notation, from now on we shall drop the subscript k — 1 in the 
notation V^fc-i for the discrete weak gradient. 
Now we introduce two forms on Vh as follows: 




s(v, w) = p ^ h T 1( \ v " Vb, W ~ W b ) d T, 

Ten 



where p > is a parameter with constant value. In practical computation, one might 
set p = 1. Denote by a s (-; •, •) a stabilization of a(-; •, •) given by 

a s (<fi;v, w) := a(<fi; v, w) + s(v, w). 



8 



Weak Galerkin Algorithm 1 . A numerical approximation for U.l\) and U.5\) 
can be obtained by seeking Uh — {u , u b } G Vh satisfying both = Q^g on dQ and 
the following equation: 

(4.7) a s (u h ;u h , v) = (/, v ), V« = {v Q , v b } G 

where Q^g is an approximation of the Dirichlet boundary value in the polynomial space 
Pk(dT n dft). For simplicity, one may take Q\>g as the standard L 2 projection of the 
boundary value g on each boundary segment. 

5. L? Projection Operators. There are two basic polynomial spaces associated 
with each element T G Th- The first one is the local finite element space Wk(T) and 
the second one is the polynomial space [Pk-i{T)] d which was utilized to define the 
discrete weak gradient V ffi in Q3.3P ; namely, the operator \? w , r ,K with r = k — 1 and 
K = T. For simplicity of discussion, we introduce the following notation 

Gfc-i(T) := [P fc _i(T)] d , 

and shall call this a local discrete gradient space. 

For each element T G Th, denote by Qq the L 2 projection from L 2 (T) onto Ph(T). 
Analogously, for each edge or flat face e G £ft, let Qb be the L 2 projection operator 
from L 2 {e) onto Pfc(e). Denote by Qh the L 2 projection onto the local discrete gradient 
space Gfc_i(T). Recall that V is the weak function space as defined by (|4.3|) . We 
define a projection operator Qh : V — > Vh as follows 

(5.1) Q h v := {Q v , QbVb}, y v = {v ,v b } eV. 

Lemma 5.1. Let Qh be the projection operator defined as in \5. Then, on each 
element T £ Th> we have 

(5.2) V m (Q^) = Q/,(V0), V^etf 1 ^). 

Proof. Using (I3.3[) . the integration by parts and the definitions of Qh and Qh, we 
have that for any r 6 Gk-i(T) 

{V w {Qh4>), t) t = -{Qo4>, V • t) t + (Qb<t>, t ■ n) 9T 
= V • r) T + {(f>, t ■ n) 9T 

= (V0, t) t = (Qfc(W), t) t 

which implies the desired relation ()5.2|) . □ 

The following lemma provides some estimate for the projection operators Qh and 
Qh- Observe that the underlying mesh Th is assumed to be sufficiently general to 
allow polygons or polyhedra. A proof of the lemma can be found in [18]. It should 
be pointed out that the proof of the lemma requires some non-trivial technical tools 
in analysis, which have also been established in |18) . 

Lemma 5.2. Let Th be a finite element partition of f2 satisfying the shape regu- 
larity assumption Al - A4. Then, for any <f> G H k+1 (n), we have 

(5.3) Yl U-Q^W 2 t+ E h 2 T \\V^-Qo^fT<Ch 2(k+ ^U\\ 2 k+1 , 
TeT h Ter h 

(5.4) - Qfc(V#)||! < Ch 2k U\\ 2 +1 . 
Ten 
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Here and in what follows of this paper, C denotes a generic constant independent of 
the meshsize h and the functions in the estimates. 

Let T be an element with e as an edge. For any function ip € iJ 1 (T), the following 
trace inequality has been proved to be valid for general meshes satisfying Al - A4 
(see [T5] for details): 



(5.5) 



y\\l<C(h^\\ip\\ 2 T + h T \\^\\ 2 T ). 



Using (|5.5p . we can obtain the following estimates. 

Lemma 5.3. Assume that l~h is shape regular. Then for any w € H k+1 (£l) and 
v = {vo,Vb} € Vh, we have 



(5.6) > ] h T l (Q Q w - Q b w, v - v b ) dT < Ch k \\w\\ k+1 \\\vl 

(5.7) J2 H Vw ~ QhVw) • n, v - v b ) dT < C7i fc ||«;|| fc+ ilHI- 
Ten 

Proof. Using the definition of Qh, (|5.5p . and (|5.3[) , we have 



^ h T 1 (Q Q w - Q b w, v -v b ) dT 



Ten 



^ h T 1 (Q w - w, v -v b )oT 



Ten 



<c(j2 ( h T 2 \\Qow - HIt + \WQoV) - w)\\ 2 T )\ ( ^ Z^^ 1 ||^o - v b \\ 2 dT \ 
\Ten ) xren ) 



<C7i fc H|fc+i|HI- 
Similarly, it follows from (|5.5I) and (|5.4I) that 



53 (a(Vw - QhVw) • n, v - v b ) dT 
Ten 

/ \ 1/2 / x 1/2 

< I 53 hrWaiVw-QhVw)^) 53 h^\ 
\Ten / VtgT^ 

<C/i*H| fc+1 ||M|. 
This completes the proof. □ 

6. On Mass Conservation. The second order elliptic equation (jl.ip can be 

rewritten in a conservative form as follows: 



V0 - Vb\\dT 



V • q = /, ? = -flVti. 

Let T be any control volume. Integrating the first equation over T yields the following 
integral form of mass conservation: 

(6.1) / q-nds= I fdT. 

JdT Jt 
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We claim that the numerical approximation from the weak Galerkin finite element 
method (|4.7[) for (| 1 . 1 [) retains the mass conservation property (16. ip with an appropri- 
ately defined numerical flux qh- To this end, for any given T £ Th, we chose in (|4.7I) 
a test function v = {v , Vb = 0} so that v — 1 on T and «o = elsewhere. It follows 
from (|4~71) that 



(6.2) / aV w Uh ■ V w vdT + ph T x / (uq — Ub)ds — / fdT. 

JT JdT Jt 

Recall that Qh is the local 1? projection onto [Pk-i(T)] d . Using the definition 
for V w v one arrives at 

/ aV w u h ■ \7 w vdT = / Qh(aV w Uh) ■ V w vdT 

JT JT 

= - I V • Q h (aV w u h )dT 

JT 

(6.3) =-/ Q h (a\7 w u h ) -nds. 

JOT 

Substituting into yields 



(6.4) / {~Q h (aV w u h )+ph- 1 (uo-u b )n}-nds= fdT, 

JdT JT 

which indicates that the weak Galerkin method conserves mass with a numerical flux 
given by 

Qh = -Qh (aV w u h ) + ph^iuo - u b )n. 

Next, we verify that the normal component of the numerical flux, namely qh ■ n, 
is continuous across the boundary of each element T. To this end, let e be an interior 
edge/face shared by two elements T\ and T 2 . Choose a test function v = {vo,Vb} so 
that vq = and v\, = everywhere except on e. It follows from (|4.7[) that 

(6.5) / aV w Uh-V w vdT—ph^ {u - u b )| Tl u b ds 

-ph?l / (u - Ub)\T 2 Vbds 

JdT 2 ne 

= 0. 

Using the definition of weak gradient p. 31) we obtain 

/ aV w ii h ■ \7 w vdT = / Qh(aV w Uh) ■ V w vdT 

JT!UT 2 JTxUT-2 

Jh(aV w u h )\ Tl ■ ni + Qh(aV w u h )\ T2 -n 2 )v b ds, 



where is the outward normal direction of Tj on the edge e. It is clear that ni +ri2 
0. Substituting the above equation into (|6.5p yields 



(-Qh(aV w u h )\ Tl + phj}(u - itf,)| Tl ni) • n x v b ds 

{-Qh{aV w u h )\ T2 + ph^{u - Ub)|r 2 n 2 ) • n 2 v b ds, 
which shows the continuity of the numerical flux qh in the normal direction. 
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7. Existence and Boundedness of WG Solutions. Let <f> £ Vh be any weak 
finite element function. A linearized version of (|4.7|) seeks Uh — {u 0l u b } e Vh 
satisfying both Ub = Q^g on dQ and the following equation: 

(7.1) a.(<f>;u h) w) = (/, «o), Vu = {t; , u 6 } £ 

It is easy to see that, for any fixed cf> € Vh, the bilinear form a s ((f>; •, •) is symmetric 
and positive definite in the weak finite element space Vh- Thus, one may introduce a 
norm in Vh as follows 

(7-2) I?; 1 := >/a s (</>;?;, v), v e V h . 

The assumptions (jl.2p and (|1 .3[) on the matrix coefficient a = a(x,r),p) imply that 
the norm ||| • |||^ are uniformly equivalent for all <fi. In particular, we shall use the norm 
arising from <f> — and denote the corresponding norm by 



The trip-bar norm ||| • ||| is an i? 1 -equivalence for finite element functions with vanish- 
ing boundary value. Moreover, the following Poincare-type inequality holds true for 
functions in V®. 

Lemma 7.1. Assume that the finite element partition Th is shape regular. Then, 
there exists a constant C independent of the meshsize h such that 

(7.3) |h||<C|H, Vt) = {vq, Vb} 6 Vjp. 

Proof. For any v — {vo, «&} 6 V®, let q 6 {H 1 (Q)] d be such that V • q = vq and 
II Q 111 < C|l u o||- To see an existence of such a function q, one may first extend vo by 
zero to a convex domain f2 which contains Q, and then consider the Poisson equation 
A "J = vo on the enlarged domain f2 and set q = V^P. The required properties of q 
follow immediately from the full regularity of the Poisson equation on convex domains. 

Recall that Qh is the L 2 projection to the space of piecewise polynomials of degree 
k - 1. Thus, 

IMI 2 - («0,V-q) T 
TeT h 

= ({uo,q-n) ar - (V«o,q)r) 
Ter h 

= X! ((«o,q-n) aT - (Vw ,Q^q)T) 
TeT h 

= V ' (^' iC 0)t ~ ( v °> i^hq) ■ n) dT + (v , q • n)ar) 

TeTh 

(7.4) = ((«o, V • (Q h q)) r - («o, (Q h q) • n) aT + (u - v b , q • n) 0T ) , 

TeTh 

where we have used the continuity of q • n across each element edge/face and the fact 
that Vb = on dil. Observe that the definition (|3.3[) of the discrete weak gradient 
implies 



(«0, V • (Q h q)) r - -(V W V, Qh^fr + (v b , (Q/,q) ■ n) aT . 
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Substituting the above identity into (|7.4p yields 

(7.5) Kf = H v ™ w > "Mr + (vo - (q - Q h q) ■ n) dT ) . 

Ten 

Using (|5.7p with w = "P, a = 1, and fc = 1, we have 



^2 (vo - Ub, (q - Q^q) • n) aT 



<d||*l| 2 |HI 



Substituting the above estimate into (|7.5I) we arrive at 

INf < IIV^H ||Q h q|| +C7i||*|| 2 |H|| 

< UVt^H ||q|| + C7»||¥|| a |H 
<C\lV w v\l |[*|| a 

< CIlV^I lluoll, 

where we have used the fact that q = and [[^Ha < C||^o|| for some constant C. 
This completes the proof of the lemma. □ 

Denote by B g .h the set of finite element functions satisfying the boundary condi- 
tion Qbg] i.e., 

Bg t h '■= {v = {vqjVi,} G Vh such that v b = Q^g on dfl}. 

The following is a result on the solution uniqueness and existence for the linearized 
problem fTTT]) . 

Lemma 7.2. The weak Galerkin finite element scheme |7.i| ] has one and only 
one solution. Moreover, there exists a constant C such that the solution of |7. 1\ ) has 
the following boundedness estimate 

(7-6) |H||| < C(||/|| + inf lll^ll). 

Proof. It suffices to show that the solution of (|7.ip is trivial if the data is homoge- 
nous; i.e., if / = g = 0. To this end, assume that the data is homogeneous. By taking 
v = Uh in (|7.ip we arrive at 



(aV w u h , V w Uh) + P X! h T X { u o ~ u b, «o - u b )or = 0, 
TeT h 



where a — a(x, <f>, V w (f>). This implies that V^it^ = on each element T and uq = u b 
on dT. It follows from V w Uh = and Q3.3P that for any q 6 [Pk-i{T)] d we have 

= {V w u h , q) T 
= -(u , V • q) T + (u b , q ■ n) dT 
= (Vu , <?)t - (uo - Ub, q ■ n) aT 
= (Vuq, q)r- 

Letting q = Vuo in the above equation yields Vito = on T € 7fe- It follows that 
Mo = const on any T £ Th- This, together with the fact that uq = u b on dT and 
Ub = on cT2, implies uq = Ub = 0. 



13 



For any tp G B g h, the difference Uh = Uh — ip is a function in satisfying the 
following equation 

a s {<j)\u h , v) = (/, u ) - a s ((f>;<ip, v), V» = {w , «&} G V,°. 

By letting v = Uh we arrive at 

\\\uh\\\l = {f, u h ) - a s {cf);ip,u h ). 

Thus, it follows from the Poincare inequality (|7.3p and the boundedness of a s (cf); •, •) 
that 

lll%lll^ < + Ml) llfihll^ 

which, together with Uh = Uh + tp and the usual triangle inequality, implies the 
designed estimate (I7.6[) . This completes the proof of the lemma. □ 

For the general nonlinear elliptic equation (|1.1[) . we have the following result on 
solution existence. 

Lemma 7.3. There exists a weak finite element function Uh £ Vh satisfying the 
weak Galerkin finite element scheme J^. 7| ). Moreover, the WG solution satisfies the 
following estimate: 

(7-7) IM, < C(||/|| + , mf mviii). 



Proof. We shall use the Leray-Schauder fixed point theorem to prove an existence 
of u h satisfying (|4.7[> . Recall that one version of the Leray-Schauder fixed point 
theorem (see for example Theorem 11.3 in [T3]) asserts that a continuous mapping F 
in K™ into itself has at least one fixed point if there exists a constant M such that 
any solution of aF(w) — w with a e [0, 1] must satisfy < M, where ||w||r™ is 

the norm of w in R™. 

For any <f> £ Vh, let £ Vh be the solution of the following linear problem: Find 
M0 = {uq, Ub} € Vh satisfying both Ub — Qb9 on dQ and the following equation: 

(7.8) a a (<ftu0, v) = (/, v ), V v = {v , v b } G V°. 

Denote by F((j>) :— the mapping from Vh into itself. It is clear that F is a 
continuous one. Assume that £/, € Vh satisfies the operator equation £/j = aF{^h) for 
some real number a G [0, 1]. This implies that £/, = crQbg on dVl and satisfies 

(7.9) 0.(&;<r- 1 & J «) - (/, v ), V v = { Vo , M G 14 . 
Multiplying both sides of (I7.9[) by er yields 

(7.10) a.(&;&, v) = (a/, u ), V v - {u , «&} G 

The estimate (|7.6p can be used to give the following estimate for the solution of (|7.10[) 

Uh% < C sup H/H + inf IIIVII). 
cre[o,i] ve-B CT9 , h 

This shows that all the conditions of the Leray-Schauder fixed point theorem are 
satisfied for the mapping F. Thus, F admits at least one fixed point Uh which is 
easily seen to be the solution of the WG finite element scheme (|4.7p . □ 
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8. Error Analysis. The goal of this section is to establish some error estimates 
for the WG finite element solution arising from (|4.7[) . Our convergence analysis 
will be established for only the linear case of (|l.ip . In other words, we shall assume 
that the coefficient matrix a = a{x, r/,p) is independent of the variables r\ and p. The 
error will be measured in two natural norms: the triple-bar norm as defined in (|7.2p 
and the standard L 2 norm. The triple bar norm is essentially a discrete H 1 norm for 
the underlying weak function. 

For simplicity of analysis, we assume that the coefficient tensor a in (jl.ll) is a 
piecewise constant matrix with respect to the finite element partition Th- The result 
can be extended to variable tensors without any difficulty, provided that the tensor a 
is piecewise sufficiently smooth. 

8.1. Error equation. Let <f> £ i? 1 (T) and v S Vh be any finite element func- 
tion. It follows from (I5.2j) . the definition of the discrete weak gradient (|3.3p . and the 
integration by parts that 

(a\7 w Q h (/), V w v) T = (aQ h (V<p), V w v) T 

= -(v , V ■ (aQ h V(j))) T + (v b , (aQ h V0) ■ n) 9T 
= (V«o, aQ h \7(j)) T - {v - v b , (aQ h Vcf>) ■ n) 9T 

(8.1) = (aVcj), Vv )t - ({aQ h V<t>) • n, v - v b ) gT . 

Testing (jl.ip by using vq of v = {vq, v b } G V® we arrive at 

( 8 - 2 ) ^2 ( aVu ' Vw o)t - ^ aVu ■ n ' v " ~ v ^ dT = 

Ten Ten 

where we have used the fact that J2reT ( a ^ u ' n > v b)dT = 0. By letting cf> = u in 
(|8.1[) . we have from combining (|8.ip and (|8.2p that 

^ {aV w Qhu, V w v)t = (/, w ) + (o(Vu - QhVu) • n, u - i%)aT- 
Adding s(QhU, v) to both sides of the above equation gives 

(8.3) a s (Q h u, v) = (/, v Q ) + ^ ( a (Vw - Q/jVu) • n, u - w b ) aT + s(Q h u, v). 

Ten 

Subtracting (|4.7p from f|8.3|) yields the following error equation 

(8.4) a a (e/,, «) = ^ ( a ( Vu _ QfcVu) ■ n, v - v b ) dT + s(Q h u, v), Vv € V,, , 

where 

e h = {e , e 6 } := {Q u - u , Q b u - u b } 

is the error between the WG finite element solution and the L 2 projection of the exact 
solution. 

8.2. Error estimates. The error equation (|8.4p can be used to derive the fol- 
lowing error estimate for the WG finite element solution. 

Theorem 8.1. Let Uh € V% be the weak Galerkin finite element solution of the 
problem il.l]) - [T75^) arising from |^.7| ). Assume that the exact solution is so regular 
that u G H k+1 (n). Then, there exists a constant C such that 

(8.5) \\u h - Q h u\\ < Ch k \\u\\ k+l . 
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Proof. By letting v — in (|8.4I) . we have 

(8.6) \\\e h \\\ 2 = ^2 ( a (Vw - QhVit) • n, e - e h ) 9T + s(Q /l u, e fe ). 

T£T h 

It then follows from (I5.6[) and (|5.7[) that 

HI ||| 2 < C/i fe ||u||fe + i|||e/ l ||| ) 
which implies (|8.5[) . This completes the proof. □ 

To obtain an error estimate in the standard L 2 norm, we consider a dual problem 
that seeks $ £ Hq(£1) satisfying 

(8.7) - V • (aV$) = e in fi. 

Assume that the usual ff 2 -regularity is satisfied for the dual problem. This means 
that there exists a constant C such that 

(8.8) ||<f||2<C||e ||. 

Theorem 8.2. Let S Vh be the weak Galerkin finite element solution of the 
problem arising from J^. 7| ). Assume that the exact solution is so regular 

that u £ H k+1 (£l). In addition, assume that the dual problem {8. 7| ) has the usual 
H 2 -regularity. Then, there exists a constant C such that 

(8.9) HQou-uoll <Ch k+1 \\u\\ k+1 . 



Proof. By testing (|8.7j) with eo we obtain 
||e || 2 = -(V-(aV$),e ) 

(8.10) = ( aV *> Ve ) T - ( aV *' n > e -e b ) dT . 

TET h TeT h 

Setting </> — $ and v — eh in (|8.1|) yields 

(8.11) (aV w Q fe $, V w e h ) T = (aV$, Ve )T - ((aQ h V$) • n, e - e & ) 0T . 
Substituting (|8~TT|) into (|8~TU|) gives 

(8.12) ||e || 2 = (aV w e h , V w Q h $) + WQfcVS - V$) • n, e - e & ) aT . 

T£T h 



It follows from the error equation (|8.4[) that 

(aV^ft, V TO <2fo$) = ^ (o(Vu - QhVu) • n, Q $ - Q b $)dT 
Ten 

(8.13) + s(g /lU , Q h $) - «(e & , Q fc $). 
By combining (|8 . 1 2[) with (|8.13p we arrive at 

(8.14) ||e || 2 = < a ( Vu - ^ Vu ) ' n > - Qb<S>)dT 

Ten 

+s(Q h u, Q h $>) - s(e h , Q h $>) 

• n, e - e 6 ) aT . 
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Let us bound the terms on the right hand side of (|8 . 14[) one by one. Using the 
Cauchy-Schwarz inequality and the definition of Qi, we obtain 



(8.15) 



E (a(V« - QhVit) • n, Q $ - Qbfyar 
Ten 

\ 1/2 

. I > \\(iiVu-Q h Vu)\\l T ) 



( E H Vl 

\TeT h 



( E IIQo*-Q6*||irj 
\TeTh / 



1/2 



<C(E IkCVu-QhVuJHlr 



1/2 



1/2 



E iiqo*-*i 

\Ten 



From the trace inequality (I5.5[) and the estimate (|5 .3[) we have 

1/2 



f E WQo$-®\\dJ) <c0\\$h 

\TET h ) 



and 



( E ii a ( Vi 



1/2 



,Vu) 



\ar 



<Ch"-?\\u\\ k+1 . 
Substituting the above two inequalities into (|8. 15|) wc obtain 



.16) 



E (o(V« - QhVu) • n, Q $ - Q b <P) dT 



Ten 



<Ch k+1 \\ u \\ k+1 m\ 



Analogously, it follows from the definition of Q^, the trace inequality (|5.5[) , and the 
estimate (15.31) that 



\s(Q h u, Q h <f>)\ < p E h T \(Qou - Qbu, Qo$ ~ Qb<$>) 



0T\ 



Ten 



<cl E ^t 1 !!^ 

\TeTh 



Ht) ( E ^IIOoS-iHlr 
/ \TeTh > 



1/2 



(8.17) <C7i K+i |MU +1 ||$|| 2 . 
The estimates ()5.6j) and (|8.5|l imply 

(8.18) |a(e h , Q h $)| < Cft||$|| 2 ||e h || < Ch k+1 ||u|| fc+1 ||*|| 2 . 
Similarly, it follows from (j5~7f and (f8~5]) that 



(8.19) 



■ n, e - e h ) 0T 



Now substituting ([STTB]) - ([87T9)l into (|8T4l) yields 

||eo|| a <C^lu||k+i||3|| a , 

which, combined with the regularity assumption (|8.8p . gives the desired optimal order 
error estimate (18.91). □ 
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9. Numerical Experiments. The goal of this section is to numerically verify 
the convergence theory for the WG finite element method (|4.7I) through some com- 
putational examples. In particular, the following issues shall be examined: 
(Nl) rate of convergence for WG solutions in various measures; 
(N2) accuracy of WG solutions on polyhedral meshes with and without hanging 
nodes. 

For simplicity, all the numerical experiments are conducted by using piecewise linear 
functions (i.e., k = 1) in the finite element space Vh as defined in (|4.5[) . 

For any given v = {vo,v b } g Vh, recall that its discrete weak gradient, V w v £ 
[Po(T)] d , is defined locally by the following equation 

(\7 w v, q) T = -K V • q) T + (v b , q ■ n) aT , Vq G [P (T)] d . 

Since q G [Po(T)] d , the above equation can be simplified as 

(9.1) (V w v,q) T = (v b ,q- n) dT . 

The error for the WG solution of (|4.7[) shall be measured in three norms defined 
as follows: 



INI 2 - £ 



I \7 w v\ 2 dx + hn 



{v - Vbfds 



or 



\v h f:= £ 
Ten 



Ml 



e£E h 



\v \ 2 dx 
\v b \ 2 ds 



(A discrete _ff 1 -norm), 
(Element-based L 2 -norm), 
(Edge-based L 2 -norm). 



9.1. Case 1: Poisson Problem on Uniform Meshes. Consider the Poisson 
problem that seeks an unknown function u — u(x, y) satisfying 

-Au = f 

in the square domain f2 = (0, l) 2 with homogeneous Dirichlet boundary condition. 
The exact solution is given by u — sin(7rx) sin(7ry), and the function / = f{x,y) is 
given to match the exact solution. 

Table 9.1 

Case 1. WG solutions and their convergence on rectangular elements. 



mcshsize h 1 


\\\Qhu - uh\\\ 


\\QhU - u h \\ 


\\Qhu - UhWe h 


4 


7.8668c-001 


1.3782e-001 


1.7244e-02 


8 


3.6731e-001 


3.5717e-002 


4.5321e-03 


16 


1.7954c-001 


9.0101e-003 


1.1362e-03 


32 


8.9221e-002 


2.2576e-003 


2.8401e-04 


64 


4.4541e-002 


5.6472e-004 


7.0995e-05 


128 


2.2262e-002 


1.4120e-004 


1.7748e-05 


0(h r ),r = 


1.0245 


1.9886 


1.9889 



Tables 19.11 and 19.21 show the rate of convergence for the corresponding WG solu- 
tions in H 1 and L 2 norms on rectangular and triangular meshes, respectively. The 
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Table 9.2 

Case 1. WG solutions and their convergence on triangular elements. 



meshsize h 1 


\\QhU - u h \\\ 


WQhU - u h \\ 


WQhU - UhWs h 


4 


1.3567c+000 


1.5399c-001 


6.5585c-02 


8 


6.8946c-001 


3.9419e-002 


1.3106c-02 


16 


3.4613c-001 


9.9131e-003 


3.0102c-03 


32 


1.7324e-001 


2.4819e-003 


7.3455e-04 


64 


8.6641e-002 


6.2072e-004 


1.8249e-04 


128 


4.3323e-002 


1.5519e-004 


4.5550e-05 


0(h r ),r = 


0.9949 


1.9925 


2.0855 



rectangular mesh is constructed by uniformly partitioning the domain into nx n sub- 
rectangles. The triangular mesh is obtained by dividing each rectangular element into 
two triangles by the diagonal line with a negative slope. The mesh size is denoted 
by h = 1/n for both the rectangular and triangular meshes. The numerical results 
indicate that the WG solution with k = 1 is convergent with rate O(h) in H 1 and 
0(h 2 ) in L 2 norms. 

9.2. Case 2: Degenerate Elliptic Problems. The second testing problem 
is defined in the square domain Q, — (0, l) 2 for the following second order partial 
differential equation 

— V • (oVti) = /, a = xy. 

Note that the coefficient a = xy > in the domain and vanishes at the origin. The 
PDE under consideration is thus elliptic, but with some degeneracy near the origin. 
The WG finite element method (|4.7|) is still applicable, and the corresponding discrete 
problem admits a unique solution. However, the convergence theory established in 
previous sections for the WG finite element method can not be applied without any 
modification. 

In our numerical tests, the exact solution is given by u = x(l — x)y(l — y), which 
corresponds to a homogeneous Dirichlet boundary condition. Like the case 1, the 
function / = f(x,y) is given to match the exact solution. 



Table 9.3 

Comparison of Convergence for Three Finite Element Schemes for a Degenerate Elliptic Problem 

WG-FEM Pi Pi DG Pi WG-FEM P P 



meshsize 
h- 1 


H -error L 2 -error 


iP-error L 2 -error 


if -error L 2 -crror 


8 


2.51c-02 1.46e-03 


3.98e-02 6.29e-03 


5.16e-02 2.01e-03 


16 


1.26e-02 3.74e-04 


3.01e-02 2.92e-03 


3.98e-02 9.30e-04 


32 


6.31e-03 9.47e-05 


2.23e-02 1.32c-03 


2.96e-02 4.02e-04 


64 


3.16e-03 2.39e-05 


1.62e-02 5.89e-04 


2.15e-02 1.70e-04 


128 


1.58e-03 6.04e-06 


1.17e-02 2.64e-04 


1.55e-02 7.17e-05 


0(h r ) 

r = 


9.97e-01 1.98e+00 


4.42c-01 1.15e+00 


4.36e-01 1.21e+00 



In Table 19.31 the column corresponding to WG-FEM Pi Pi refers to the com- 
putational results obtained from the numerical scheme (I4.7|) with piecewise linear 
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functions on each element and its edges. The column corresponding to DG Pi is the 
result arising from the interior penalty method with piecewise linear functions. The 
last column corresponds to results from the weak Galerkin method detailed in [17) 
with piecewise constants. These three methods were chosen for comparison because 
they have the same rate of convergence in theory when the error is measured between 
the finite element solution and a certain interpolation of the exact solution. 

The computational results indicate that the new WG-FEM scheme (|4.7|) presented 
and analyzed in the present paper has optimal order of convergence in both H 1 and 
L 2 , while the other two converges with significantly lower orders. The H 1 norm in 
the table refers to discrete equivalence for each respective scheme. 

9.3. Case 3: WG-FEM on Deformed Rectangular Meshes. We solve 
the same problem as in Case 1 on deformed rectangular meshes. We start with 
an initial deformed rectangular mesh, shown as in Figure |9~T1 (Left). The mesh is 
then successively refined by connecting the baryccntcr of each (coarse) element with 
the middle points of its edges, as shown in the dotted line in Figure |9~T1 (Right). 
The numerical results are presented in Table 19.41 which show an optimal order of 
convergence in various norms. 




0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 



Fig. 9.1. Case 3: An initial mesh (Left) and its refinement (Right). 



Table 9.4 

Case 3. Convergence rate on deformed rectangles. 



meshsize h 


\\\QhU - u h \\\ 


WQhU - UhW 


\\QhU - UhWs h 


2.8790e-01 


2.3056c+00 


3.0235e-01 


8.2633e-02 


1.4395e-01 


1.1673e+00 


7.8108e-02 


2.1396c-02 


7.1974e-02 


5.8473c-01 


1.9652c-02 


5.3912e-03 


3.5987e-02 


2.9241c-01 


4.9203c-03 


1.3503e-03 


1.7993e-02 


1.4619e-01 


1.2445e-03 


3.3774e-04 


8.9967e-03 


7.3095e-02 


3.1112e-04 


8.4445e-05 


0{h r ),r = 


0.9828 


1.9618 


1.9893 



9.4. Case 4: WG-FEM on Meshes with Hanging Nodes. We solve the 
same problem as in Case 1 on deformed rectangular meshes with hanging nodes in 
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the finite element partition. The initial mesh is shown as in Figure (Left). The 
mesh on the right in Figure HOI is generated by following the same uniform refinement 
procedure as described in Case 3. It should be pointed out that the initial mesh has 
a hanging node in the usual definition. 

For the finite element partition Th with hanging nodes, the WG finite element 
method (|4.7j) must be modified as follows. For edge containing hanging nodes, the 
edge shall be further partitioned into smaller segments by using the hanging nodes. 
Then the corresponding finite element space defined on this edge will be piecewise 
linear functions with respect to the new partition; the finite element space on each 
element remains unchanged. For example, in Figure [931 the elements K\, K2, and K3 
share one hanging node M. In the WG finite element method, the edge AB needs to 
be divided into two pieces: AM and MB. The corresponding finite element function 
vi, is taken as a piecewise linear function on edges AM and MB. We point out 
that the refinement method adopted here may produce elements around the hanging 
node which are not shape regular as defined in Section 3. The numerical results are 
presented in Table l9~5l Readers are encouraged to draw conclusions from this table. 

Table 9.5 

Case 4- WG solutions and their convergence on deformed rectangular elements with hanging 
nodes. 



meshsize h 


\\\QhU - u h \\\ 


WQhu - u h \\ 


WQhu - Uh\\e h 


4.2512c-01 


3.8064c+00 


9.1677e-01 


2.5509e-01 


2.1256c-01 


2.2593c+00 


3.2970e-01 


7.3355e-02 


1.0628e-01 


1.2308e+00 


9.2302e-02 


2.2037e-02 


5.3140e-02 


6.3470e-01 


2.3300e-02 


6.1413e-03 


2.6570e-02 


3.2104e-01 


5.8376c-03 


1.7649e-03 


1.3285e-02 


1.6129c-01 


1.7094c-03 


5.1822c-04 


0(h r ),r = 


0.9201 


1.8508 


1.7912 




0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 



Fig. 9.2. Case 4: Mesh level 1 (Left) and Mesh level 2(Right). 
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